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Abstract: The switch-Uke character of gene regulation has motivated the use of hybrid, 
discrete-continuous models of genetic regulatory networks. While powerful techniques for 
the analysis, verification, and control of hybrid systems have been developed, the specifici- 
ties of the biological appHcation domain pose a number of challenges, notably the absence 
of quantitative information on parameter values and the size and complexity of networks 
of biological interest. We introduce a method for the analysis of reachability properties of 
genetic regulatory networks that is based on a class of discontinuous piecewise-afHne (PA) 
differential equations well- adapted to the above constraints. More specifically, we introduce 
a hyperrectangular partition of the state space that forms the basis for a discrete abstraction 
preserving the sign of the derivatives of the state variables. The resulting discrete transition 
system provides a conservative approximation of the qualitative dynamics of the network 
and can be efficiently computed in a symbolic manner from inequality constraints on the 
parameters. The method has been implemented in the computer tool Genetic Network 
Analyzer (GNA) , which has been appHed to the analysis of a regulatory system whose func- 
tioning is not well-understood by biologists, the nutritional stress response in the bacterium 
Escherichia coli. 

Key- words: Piecewise-afRne differential equations, qualitative analysis, discrete abstrac- 
tion, hybrid systems, genetic regulatory networks, systems biology, nutritional stress re- 
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Analyse symbolique d'atteignabilite de reseaux de 
regulation genique par abstractions qualitatives 

Resume : Le fonctionnement binaire allume/ eteint des genes a motive rutilisation de 
modeles hybrides pour representer les reseaux de regulations geniques. Bien que des techni- 
ques efBcaces aient ete developpees pour I'analyse, la verification et le controle des systemes 
hybrides, les specificites des systemes biologiques consideres, en particulier I'absence d'infor- 
mations quantitatives sur les valeurs des parametres et la taille et la complexite des reseaux 
d'intcret biologique, posent un certain nombre de problemes nouveaux. Nous proposons une 
methode pour I'analyse des proprietes d'atteignabilite des reseaux de regulations geniques 
basee sur une classe de modeles bien adaptes aux contraintes mentionnees ci-dessus, uti- 
lisant des equations difFerentielles discontinues et afRnes par morceaux. Plus precisement, 
nous introduisons une partition hyperrectangulaire de I'espace d'etat, servant a definir une 
abstraction discrete qui preserve les signes des derivees des variables d'etat. Le systeme de 
transitions discret ainsi defini est une approximation conservative de la djmamique du reseau 
et pent etre calcule symboliquement de fagon efiicace a partir de contraintes d'inegalite sur 
les parametres. Cette methode a ete implementee dans I'outil Genetic Network Analyzer 
(GNA), et utilisee pour ranalyse d'un systeme biologique encore mal compris, celui de la 
rcponsc au stress nutritionncl chcz la bactcric Escherichia coli. 

Mots-cles : Equations differentielles afRnes par morceaux, analyse qualitative, abstraction 
discrete, systemes hybrides, reseaux de regulations geniques, biologic des systemes, reponse 
au stress nutritionnel, Escherichia coli 
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1 Introduction 

The functioning and development of living organisms is controlled on the molecular level 
by networks of genes, proteins, small molecules, and their mutual interactions, so-called 
genetic regulatory networks. The dynamics of these networks is hybrid in nature, in the 
sense that the continuous evolution of the concentration of proteins and other molecules is 
punctuated by discrete changes in the activity of genes coding for the proteins. The switch- 
like character of the dynamics of genetic regulatory networks has attracted much attention 
from researchers on hybrid systems {e.g., dEl, 21, H, 36|). 



While powerful techniques for the analysis, verification, and control of hybrid systems 
have been developed (see 0, [iH for reviews), the specificities of the biological application 
domain pose a number of challenges [lil . l4a |. First, most genetic regulatory networks of 
interest consist of a large number of genes that are involved in complex, interlocking feedback 
loops. Second, the data available on both the structure and the dynamics of the networks is 
currently essentially qualitative in nature, meaning that numerical values for concentration 
variables and kinetic parameters describing the interactions are generally absent. The above 
characteristics require hybrid-system methods and tools to be upscalable and capable of 
dealing with qualitative information. 

In this paper, we will show that a class of piecewise-affine differential equation (FADE) 
models, originally introduced by Glass and Kauffman in the seventies [3o|, is particularly 
suitable for deahng with the above challenges. The properties of these FADE models have 
been well-studied in mathematical biology [H, H, 0, IH, S 113, SI, H SI. 



The variables in the FADE models are the concentrations of the proteins encoded by the 
genes, while step functions account for the discrete changes in gene activity occasioned by 
the regulatory interactions. On a formal level, the FADE models are related to a class of 
asynchronous logical models proposed by Thomas and colleagues [i^, fiS]. FADE models 
and their logical relatives have been used for the study of a number of prokaryotic and 
eukaryotic regulatory networks (see [i^ for a review and references). 

The particular form of the FADE models allows the continuous state space to be par- 
titioned into hyperrectangular regions with equivalent qualitative dynamics, preserving the 
sign pattern of the time derivatives of the solutions. We exploit this partition by associ- 
ating to each FADE model a continuous transition system having equivalent reachability 
properties. By means of a discrete or qualitative abstraction 0, UQ, MS, Si, the con- 
tinuous transition system is turned into a discrete transition system providing a compact 
and qualitative description of the dynamics of the continuous system. Formally, there exists 
a simulation relation between the continuous and discrete transition systems, that is, the 
latter is a conservative approximation of the former. We show that the discrete transition 
system is invariant over large ranges of parameter values and can be computed in a symbolic 
manner from inequality constraints. Moreover, it is possible to design tailored algorithms 
for computing the discrete transition system, which scale up to large and complex genetic 
regulatory networks. 

The paper continues our previous work on the qualitative analysis of the dynamics of 



genetic regulatory networks by means of FADE models |2ll . l22l | . The main novelty of this 



RR n° 0123456789 



4 



G. Batt et al. 



contribution, a preliminary version of which was presented in |1Q|, is the use of a more fine- 
grained discrete abstraction, preserving the derivative sign pattern. This is an important 
feature for the experimental vaHdation of models of genetic regulatory networks, since mea- 
surements of gene expression often result in observations of changes in the sign of derivatives. 
Notwithstanding this increase in precision, the refinement does not threaten the applicabil- 
ity of the approach to large and complex networks. This is illustrated by the analysis of a 
network that is not yet well understood by biologists, the network controUing the carbon 
starvation response of the bacterium Escherichia coli. The appHcation of the method has 
led to novel insights into how the adaptation of cell growth to nutritional stress emerges 
from the network of molecular interactions. 

The paper is organized as follows. In Sections [2] and [3] we specify the PADE models 
of genetic regulatory networks in detail and discuss their mathematical properties, paying 
special attention to compHcations arising from the discontinuities in the righthand-side of the 
differential equations. In Section |4] we define a qualitative abstraction of the dynamics of PA 
systems, based on a hyperrectangular partition of the state space. Section [5] introduces rules 
to actually compute the discrete transition system induced by the qualitative abstraction 
and the implementation of the rules in a computer tool called GNA. In Section[6]we illustrate 
the application of the method to the analysis of the E. coli network. In the final section we 
present our conclusions and discuss the results in the context of related work. 



2 PADE models of genetic regulatory networks 

The dynamics of genetic regulatory networks can be described by a class of piecewise-affine 
differential equations (PADE) models of the following general form 3^, i3| : 

X = h{x) ^ f{x) ~ g{x)x, xeVt\Q, (1) 

where x = {xi, . . . ,Xn)' is a vector of cellular protein concentrations, / — (/i, • ■ • , /«)', 
g — diag((7i, . . . , .g„), Q. C M"q is a bounded n-dimensional state space region, and 8 a 
zero-measure subset of Vl (see below). The rate of change of each protein concentration Xi, 
i G {1, . . . , rt}, is thus defined as the difference of the rate of synthesis fi{x) and the rate of 
degradation gi{x) Xi of the protein. 

The function : 0\8 — > R>o expresses how the rate of synthesis of the protein encoded 
by gene i depends on the concentrations x of the proteins in the cell. More specifically, the 
function fi is defined as 

Mx) = Y.K\biix), (2) 

where k\ > is a rate parameter, b\ : il\Q {0, 1} a piecewise-constant regulation function, 
and Li a possibly empty set of indices of regulation functions. The function gi expresses the 
regulation of protein degradation. It is defined analogously to fi, except that we demand 
that gi is strictly positive. In addition, in order to formally distinguish degradation rate 
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parameters from synthesis rate parameters, we will denote the former by 7 instead of k. 
Notice that with the above definitions, ft. is a piecewise-affine (PA) vector- valued function. 

A regulation function 6- describes the conditions under which the protein encoded by 
gene i is synthesized (degraded) at a rate (7' Xi). It is defined in terms of step functions 
and is the arithmetic equivalent of a Boolean function expressing the logic of gene regulation 
(sol . I47I . More precisely, the conditions for synthesis and degradation are expressed using 
the step functions s+,s~: 



1, if x-j > 
0, if Xj < 



'ixj,9j) = 1 - s+{xj,9j), 



(3) 



where Xj is an element of the state vector x and 6j a constant denoting a threshold concen- 
tration. Notice that the step functions are not defined at the thresholds. 

Figure [Ija) gives an example of a simple genetic regulatory network consisting of two 
genes, a and 6. When a gene (a or 6) is expressed, the corresponding protein (A or B) 
is synthesized at a specified rate (kq or Hh). Proteins A and B regulate the expression of 
genes a and b. More specifically, protein B inhibits the expression of gene a, above a certain 
threshold concentration 0b, while protein A inhibits the expression of gene 6 above a thresh- 
old concentration 0}^, and the expression of its own gene above a second, higher threshold 
concentration 0^. The degradation of the proteins is not regulated and proportional to the 
concentration of the proteins (with degradation parameters -fa or 76). The FADE model of 
this network is shown in Figure [Ufb) . 



B ^ 



(a) 



Xa = KaS {Xa,0l)s {xf,, 0b) - Xa, 
Xb KbS''{Xa,0l) - JbXb- 

(b) 



Figure 1: (a) Example of a genetic regulatory network of two genes (a and 6), each coding for 
a regulatory protein (A and B). See Figure[6]for the legend, (b) FADE model corresponding 
to the network in (a). 



The use of step functions s^{xj, 0j) in |(T]) gives rise to mathematical complications, be- 
cause the step functions are undefined and discontinuous at Xj — 0j. Therefore, h is unde- 
fined and may be discontinuous on the threshold hyperplanes 9 = lJie{i n} heli Pi}^^ ^ 



n 



6*/}, where pi denotes the number of thresholds for the protein encoded by gene 
i. In order to deal with this problem, we can follow an approach originally proposed by 
Filippov ^!7\ and widely used in control theory. It consists in extending the differential 
equation x — h{x), x G \ 0, to the differential inclusion 



X £ K{x), with K{x) — co({ 



lim /i(y)}), x(^n, 



(4) 
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where W{P) denotes the smallest closed convex set containing the set P, and 
{\imy^x,y^e h{y)} the set of all Hmit values of h{y), for y ^ <d and y x. This ap- 
proach has been applied in the context of genetic regulatory network modeling by Gouze 
and Sari [|32|. 

In practice, K{x) may be difficult to compute because the smallest closed convex set 
can be a complex polyhedron in fl. We therefore employ an alternative extension of the 
differential equation: 

X e H{x), with H{x) = rect{{ lim h{y)}), x e il, (5) 

where rect[P) denotes the smallest closed hyperrectangular set containing the set P. The 
advantage of using red is that we can rewrite H{x) as a system of differential inclusions 
Xi € Hi{x), i G {1, . . . , n}. Notice that H{x) is an overapproximation of K{x), for all x e f2. 

Formally, we define the PA system S as the triple {fl,Q,H), that is, the set-valued 
function H given by and defined on the n-dimensional state space with 9 the union 
of the threshold hyperplanes. A solution of the PA system E on a time interval / is a 
solution of the differential inclusion ^ on /, that is, an absolutely-continuous vector- valued 
function £,{t) such that ^(t) € H{^{t)) almost everywhere on /. In particular, the derivative 
of ^{t) may not exist, and therefore £,{t) £ H{(_{t)) may not hold, if ^ reaches or leaves 9 at 
t. 

For all xq G and t G M>o U {oo}, Ss(a;o,T) denotes the set of solutions £,{t) of 
the PA system S, for the initial condition ^(0) — xo, and t S [0,r]. The existence of at 
least one solution ^ on some time interval [0, r], r > 0, with initial condition ^(0) = xq is 
guaranteed for all xq in However, there is, in general, not a unique solution. The 

set = UaoGSi T>o ^^i^Oi t) is the set of all solutions, on a finite or infinite time interval, 
of the PA system S. We restrict our analysis to the solutions in Ss that reach and leave a 
threshold hyperplane finitely-many times. The dynamics of S is thus defined by the set of 
solutions Ss. 

3 Mathematical analysis of PA systems 
3.1 Mode domains 

The dynamical properties of the solutions of S can be analyzed in the n-dimensional state 
space hyperrectangle 51 = ili x . . . x f7„, where fli = {xi e M | < < maxi} and maxi 
denotes a maximum concentration for each protein, i € {l,...,n}. In particular, we set 
maXi > meLX^fzQ\Q fi{x)/ gi{x). It is not difficult to show that under this condition 17 is a 
positively invariant set. 

For subsequent use, we introduce the notion of hyperrectangular partition of a set i? = 
i?i X . . . X Rn C SI. The partition is induced by sets of hyperplanes orthogonal to one 
of the axes Xi, i £ {1, . . . ,rt}. More precisely, the hyperplanes orthogonal to the x^-axis 
are given by the finite sets A^ C fli, where every X € Ai corresponds to a hyperplane 
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{x ^ \ Xi = A}. The hyperrectangular partition of R induced by A = {Ai, . . . , A„} is 
defined asV = PiX . . .xV^, where Pi, i & {1, ... ,n}, is the interval partition of Ri induced 
by Ai. That is, Pi is the partition of minimal cardinaUty of Ri, such that for every P GPi, 
either P = {A} for some A S A.^, or P is an interval containing no A G Aj. As an example, 
consider the interval Ri = {ri,r2] and A^ = {Ai, A2}, with ri < Ai < r2 < A2. The interval 
partition of i?i induced by A.^ equals Pi = {(ri, Ai), {Ai}, (Ai, r2]}. 

The {n — l)-dimensional threshold hyperplanes introduced in Section [2] define a hyper- 
rectangular partition of fi. 

Definition 1 (Mode domain partition). M is the hyperrectangular partition of induced 
by {9l, . . . ,6f'}, i G {!,...,«}. The sets M G M are called mode domains. 

The reason for speaking of mode domains will become clearer below, when we show that 
the regulation of gene expression is identical in every M G , and thus corresponds to a 
regulatory mode of the system. 



maxf, 
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Figure 2: (a) Mode domain partition of the state space for the model of Figure Hfb). 
(b) Focal sets and vector fields associated with the mode domains to , and M^^. 



Figure [D^a) shows the mode domain partition of the two-dimensional state space of the 
example network. We distinguish between mode domains like and M*", which are located 
on (intersections of) threshold hyperplanes, and mode domains like , which are not. The 
former domains are called singular mode domains and the latter regular mode domains. We 
denote by Mr and Ms the sets of regular and singular mode domains, respectively. Note 
that M = Mr^Ms. 

We next introduce some basic topological concepts. For every hyperrectangular set, R C 
ri, of dimension k, k & {0, . . . , n}, we define the supporting hyperplane of R, supp{R) C Q., 
as the fc-dimensional hyperplane containing R. The boundary of R in supp{R), denoted 

o o 

by dR, is defined as the set R\R, where R is the closure of R and R its interior, both 
relative to supp[R). Suppose that M is a singular mode domain, i.e. M G Ms. Then 
R{M) is defined as the set of regular mode domains M' having M in their boundary, i.e. 

R{M) = {AP G Mr I M C dhP}. 
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Using the definition of the differential inclusion |[5|), it can be easily shown that in a 
regular mode domain M, 

H{x) = {^^^ - ly^^ x}, X e Af, (6) 
where /i*-'^ is a vector of (sums of) positive synthesis rate constants and z/*^ = diag(i^f-'^, . . . , u^') 



a diagonal matrix of (sums of) positive degradation rate constants (22| ■ We define the focal 
set 

■^{M) = {V(M)}, (7) 
where %Ij(M) — {v^'^)^^ /i*^, and repeat the following standard result (3o| . 

Lemma 1. Let M G Air- Every solution ^{t) G Ss in M monotonically converges towards 
^'(M). 

Proof. From |(6]) it follows that the solutions in a regular mode domain M are given by 
^{t) — ip{M) ~ *{ip{M) — xo), for initial conditions xo £ M. As a consequence, ^{t) 
monotonically converges towards vE'(Af) = {%1){M)} as long as ^(t) e M . □ 

We will make the generic assumption that the focal sets '^{M), for all AI G M.r, are 
not located in the threshold hyperplanes 9. Figure [2ljb) shows the focal sets of four regular 
mode domains (M\ M^, and M"). In the case of M", we see that ^'(M") C M", 
so that xl)[M^^) is an asymptotically stable equihbrium point of S. 

In a singular mode domain, the right-hand side of the differential inclusion |[5|) reduces 
to the following set: 

H{x) = TEEiiifj.'^'' ~ v^^' X I M' e i?(M)}), X e M. (8) 

The focal set associated with the domain is now defined as 

*(M) = supp{M) n 7irf({V'(M') I M' e i?(M)}), (9) 

which is generally not a single point in higher-dimensional systems. As will be shown below, 
two different cases can be distinguished. If ^'(M) is empty, then every solution passes 
through M instantaneously (Lemma [2]). If not, then some (but not necessarily all) solutions 
arriving at M will remain in M for some time, sliding along the threshold hyperplanes 
containing M (Lemma[3]). In the latter case, weaker monotonicity properties hold (Lemmas[3] 
andS]). 

Lemma 2. Let M e Ms- Every solution ^(t) e Ss arriving at M instantaneously crosses 
the domain if and only if ^'(M) = {}. 

Proof. We first prove sufficiency. From the definition of ^'(Af) in a singular mode do- 
main and '^{M) = {}, it follows that for some i G {!,..., n} such that Mi is included 
in a threshold hyperplane, the intersection of suppi{M) = {0^'}, h G and 
rect{{'4)i{M') I M' G i?(A/)}) is empty. As a consequence, either Vi(^') > 6*!" for all 
M' e R{M), or ip^{M') < 6'-* for all M' £ R{M). From ^ and the hyperrectangular shape 
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of H{x), we obtain H,{x) = rect{{vf (MM') ~ e*-*) | M' e R{M)}), for all x e M. As 
a consequence, either mini?i(x) > or maxiJi(x) < 0, so that ^ and no solution 

with ^i(t) = exists in M. Therefore, every solution arriving at M instantaneously crosses 
the domain. 

Necessity is proven by contraposition, by an argument similar to that given for sufficiency. 
From \I>(M) ^ {} we show that G Hi{x), for every x £ M and i € {1, . . . ,n} such that 
Mi is included in a threshold hyperplane. We can then define solutions S,{t) that satisfy 
^i{t) = on some time interval, that is, solutions that do not instantaneously cross M. □ 

Lemma 3. Let M £ Ms and '^{M) ^ {}. For every solution S^{t) E Ss in M, and every 
i G {1, . . . ,n} such that Mi is included in a threshold hyperplane, it holds that ^i{t) — 0. 
For all other i, £,i{t) monotonically converges towards ^'i(M). 

Proof. In order to remain in M, ^(i) must slide along the threshold hyperplanes containing 
M. Therefore, £,i{t) = for every x G M and « G {1, . . . , n} such that Mi is included in a 
threshold hyperplane. For all other i, we must have ^i{t) G Hi{^{t)) almost everywhere, by 
the definition of solutions of PA systems. Moreover, Hi{£^{t)) = rect{{vf^ {tpi{M') — £,i{t)) \ 
M' G RiM)}). Let Ciit) ^ «'i(M), and ^^{t) < ip^iM') for all M' G i?(M) (the case 
^i{t) > ipi{M') goes analogously). It follows that min Hi{^{t)) > 0, and hence that ^i{t) > 0, 
provided that ^(t) G H{£^(t)). We thus infer that ^i{t) monotonically converges towards 
*»(M). □ 

Lemma 4. Let M G Ms and \f(M) 7^ {}. For every tp G *(M) and a; G M, there exists a 
solution £^{t) G Ss in M, with ^(0) = x, such that ^{t) monotonically converges towards V'- 

Proof Let M+(i) = argmaxM'efi(M) tpi{M') and M~(z) = argminM'ei?(Af) A{M'), for 
all i G {1, • • ■ Notice that due to the hyperrectangular shape of ^(M), we can write 
ipi = Oii 'ipi{M~{i)) + (1 — Ui) 'ipi{M~^ (i)), for some Ui G [0,1]. Moreover, after defining 

Pi G [0, 1] such that Pi — ai vf^ '■*V(Q^i + (1 — a^) vf^ ^^^), we introduce /i^ and Vi 

such that 

That is, ipi is written as the quotient of /i^ and fi, which in turn can be defined as a convex 
combination of /if^ , and t'f ^ , lyf' , respectively. 

Now, construct a function ^(t) = ^ — e^"^* (ip — x). This yields ^i(i) = Vi e^"'* (V'i — a^i) = 
J^i (^i — ^i(i)) = l^i — Vi £,i{t), for every i G {1, . . . , n}. Following the definition of /i^ and 
in (fTOl) . this gives 

= A (Mf"^^^ - -f"^'^6(0) + (1 - A) (/.f - -f'<^^e.(i)). (11) 

For ^(i) to be a solution, it must satisfy ^i(t) G Hi{^{t)) almost everywhere, for every 
i G {1, . . . , rt}. It can be directly verified that the expression for £_i{t) in with Pi G [0, 1], 
is included in Hi{£_{t)) as defined by iJH). Hence, the condition is satisfied and ^(t) is a solution 
of S in M. Moreover, ,^(0) — x and S,{t) monotonically converges towards ip. □ 
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In the sequel, domains which every solution crosses instantaneously will be called in- 
stantaneous domains, whereas domains in which at least some solutions remain for some 
time will be referred to as persistent. Figure [^b) shows two singular mode domains, 
and M'^. is an instantaneous mode domain (^'(M^) — 0), whereas is a persistent 
mode domain (^(M"*) — {{0^,0)'}) in which solutions slide along the threshold plane. In 
this simple example, it is intuitively clear how to define the fiow in M'^, given the dynamics 
in and M^. The use of differential inclusions as described above makes it possible to 
define the fiow in singular domains in a systematic and mathematically proper way. 

The fact that every mode domain is associated with a unique focal set has provided 



the basis for a discrete abstraction criterion employed in our previous work |15l . l21l . l22l |. 
However, this criterion disregards the fact that the system does not always exhibit the same 
qualitative dynamics in different parts of a mode domain, in the sense that the sign pattern 
of the derivatives of the solutions may not be unique. Consider a solution ^{t) e in M^^ 
in Figure[2llb): depending on whether ^t{t) is larger than, equal to, or smaller than the focal 
concentration Kb/jb in M^^, ^b will be decreasing, steady, or increasing. As a consequence, 
if we abstract the domain M^^ away into a single discrete state, we will not be able to 
unambiguously infer that solutions entering this domain from are increasing in the Xb- 
dimension. This may lead to problems when comparing predictions from the model with gene 
expression data, for instance the observed variation of the sign of ib- Today's experimental 
techniques, such as quantitative RT-PCR, reporter genes, and DNA microarrays, usually 
produce information on changes in the sign of the derivatives of the concentration variables. 



3.2 Flow domains 

The mismatch between the description levels of the mathematical analysis and the exper- 
imental data calls for a finer partitioning of the state space, which can then provide the 
basis for a more adequate abstraction criterion. Along these lines, the regular and singular 
mode domains distinguished above are repartitioned by means of the {n — l)-dimensional 
hyperplanes corresponding to the focal concentrations. 

Definition 2 (Flow domain partition). is the hyperrectangular partition of a mode 
domain M eM induced by {V',;(M)}, if M is regular, and by {'^p^iM') \ M' e R{M)}, if M 
is singular, i Q {1, . . . ,n}. T> — UMeA^I?*^ is a partition of and the sets D E V are called 
flow domains. 

The reason for speaking about fiow domains is that, as will be seen below, in every 
D E V the fiow of the system is qualitatively identical. The partitioning of the state space 
into 27 fiow domains is illustrated for the example system in FigurelHKa). Every fiow domain 
is included in a single mode domain, a relation captured by the surjective function mode: 
V Ai, defined as mode[D) = M, if and only if D C M. Similarly, the function flow: 
n ^ V denotes the surjective mapping that associates a point in the state space to its fiow 
domain: flow{x) — D, if and only if x € D. 

The repartitioning of mode domain M^^ leads to six fiow domains (Figure [Slja)). The 
finer partition guarantees that in every fiow domain of M^^, the derivatives have a unique 
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Figure 3: (a) Flow domain partition of the state space of the model in Figure [ijb). (b) State 
transition graph of the corresponding qualitative transition system. For the sake of clarity, 
self-transitions are represented by dots and transition labels are omitted, (c) Parameter 
inequality constraints for which the graph in (b) is invariant. 



sign pattern. In D^^-^, for instance, the Xa-derivative is negative and the a;f,-derivative is 
positive, whereas in D^^-^ both derivatives equal zero (in fact, D^^-^ coincides with ip{M^^) 
and is an equilibrium point of the system). This property is true more generally. Consider a 
point a; in a flow domain D eV. We denote by S{x) £ 2'^^^^°'^^" the derivative sign pattern 
at X, that is, the set of derivative sign vectors of the solutions in D passing through x. More 
formally, 

S{x) = {signiiir)) \ ^{t) G 5s in D, and 3t > : ^(t) = x and C(t) G H{^{t))}. (12) 
This gives the following central result. 

Theorem 1 (Qualitatively-identical dynamics). For all x,y & D and D &T>, S{x) — S{y). 

Proof. Note that due to the hyperrectangular shape of H{x), S{x) can be decomposed into 
5*1(2;) X ... X Snix), where Si{x) denotes the set of the ith components of the derivative 
sign vectors of the solutions in D passing through x, i G {1, . . . ,n} (idem S{y)). Define 
M = mode(D) and distinguish the cases (a)-(c). 

(a) M G Mr- Suppose 1 G Si{x), but 1 ^ Si{y), for any i G {1, . . . ,n} (for and -1 the 
argument is similar). By definition of Si{x) this means that ^^(t) > for some S,{t) G in 
D and r > 0, such that ^(r) = x. On the other hand, if>i{cj) < for all ip{t) G in D and 
cr > 0, such that ip{a) = y. Given that ^^(t) G Hi{^{T)) = {vf^ (AiM) ~ Cj(t))}, due to 
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the definition of S{x), it follows that ^{(t) = Xi < ipi{M). Similarly, (pi{(j) = yi > %lJi{M). 
But then, by Definition [11 x and y are not located in the same flow domain D, contrary to 
the assertion of the theorem. Therefore, 1 G Si{x) impHes 1 G Si{y). The converse is shown 
in the same way. 

(b) M G AAs and ^'(Af) — {}. No solution remains in D, so by definition >5'(a;) = 

S{y) - {}. 

(c) M G Ms and '^{M) ^ {}. For every i G such that A is located in 
a threshold or focal hyperplane, we must have ^i(i) = for solutions remaining in D. 
Consequently, Si{x) — {0} (idem Si{y) = {0}). For all other i, suppose 1 G Si{x), but 
1 ^ Si{y) (the argument for and —1 is similar). It follows that ^i{T) > for some £,{t) G 

in D and t > 0, such that ^(r) = x. However, (/3i((T) < for all (p{t) G in D and u > 0, 
such that ifia) = y. FVom iir) G ^^^(^(t)) = (V'»(M') - ^.(t)) | M' G i?(M)}), 

we conclude that ^iir) = Xi < 2pi{M'), for some M' G R{M). Similarly, we must have 
(^.,(cr) G H,{ip{a)) = (V'i(M') - (^,(cr)) | Af G i?(Af)}). Now, if there were 

some Af G R{M) such that </5i(cr) = yi < ipi{M'), then by Lemma |4] there would exist 
a solution such that </ji(cr) > 0. Since this contradicts the assumption, we conclude that 
(pi{cr) — yi > ipi{M'), for all Af G R{M). This implies, again by Definition [2l that x and 
y are not located in the same flow domain D. Therefore, 1 G Si{x) impHes 1 G Si{y) (and 
conversely) . □ 

Notice that the deflnition of S" as a set is a direct consequence of the use of differential 
inclusions. Since the solutions of differential inclusions are not unique, several solutions 
may pass through x and their derivatives may have a different sign in some dimension 
i G {1, . . . , rt}. This situation does not occur in our two-gene example. 

Theorem [T] suggests that the partition of the state space introduced in this section can be 
used as an abstraction criterion better-adapted to the available experimental data on gene 
expression. This idea will be further developed in the next section. 

4 Qualitative abstraction of the dynamics of PA systems 
4.1 Qualitative PA transition systems 

As a preparatory step, we deflne a continuous transition system having the same reachability 
properties as the original PA system E. Consider x € D and x' G D' , where D,D' G V 
are flow domains. If there exists a solution ^(t) of S passing through x at time t G ]R>o 
and reaching x' at time r' G M>o U {oo}, without leaving DUD' in the time interval [r, r'], 
then the absolute continuity of ^{t) implies that D and D' are either equal or contiguous. 
More precisely, one of the following three cases holds: D = D' , D C dD' , or D' C dD. 
We consequently distinguish three corresponding types of continuous transitions: internal, 

denoted by x x', dimension increasing, denoted by x '^-^ x' , and dimension decreasing, 

denoted by x '''^—^ x' . The latter two terms refer to the increase or decrease in dimension 
of the flow domain when going from D to D' . This leads to the following deflnition: 
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Definition 3 (PA transition system). S-TS = (fi. L, H, — >, |=) is the transition system 
corresponding to the PA system S = {ft, 9, H), where: 

1. is the state space; 

2. L — {int,dim^ ,dim^} is a set of labels denoting the three different types of transi- 
tions; 

3. n = {Dstgn = S \ S e 2^-'^ °''^^"} is a set of propositions describing the derivative 
sign pattern of the concentration variables; 

4. is the transition relation describing the continuous evolution of the system, defined 

hj n X L X il, such that x ^ x' if and only if there exists £,{t) € Ss and r, r', 
such that < T < t', ^(t) = x, £,{t') — x' , and 

• if ^ = int, then for all t £ [t, t']: ^{t) G flow{x) = flow{x'), 

• ii I — dim~^ , then for all t G (t, t']: £^(t) G flow{x') ^ flow{x), 

• ii I = dim~ , then for all i G [t, r'): ^(i) G flow{x) / flow{x'); 

5. \= is the satisfaction relation of the propositions in 11, defined by O x 11, such that 
X 1= Dsign = S" if and only if = S{x). 

The satisfaction relation \= thus associates to each point x in the state space a qualita- 
tive description of the dynamics of the system at x. We define any sequence of points in 
rj, . . . , x"*), m > 0, as a run of S-TS if for all j G {0, . . . , to — 1}, there exists some 

/ G L such that x^ It is not difficult to show that a PA system E and its corre- 

sponding PA transition system S-TS have equivalent reachability properties (see the proof 
in Appendix [Aj . 

Proposition 1 (Reachability equivalence). For all x,x' G il, there exists a solution ^(i) of 
E and t, t', such that < r < r', ^(r) = a;, and ^(r') = x' if and only if there exists a run 
(x°, . . . , x™) of E-TS such that x^ = x and = x'. 

The continuous PA transition system has an infinite number of states and transitions, 
which makes it difficult to verify dynamical properties of interest, for instance by means 



of conventional tools for model checking [l7|. However, we can define a discrete transition 
system, with a finite number of states and transitions, which preserves important properties 
of the qualitative dynamics of the system. In order to achieve this, we introduce the equiva- 
lence relation ~a CQ,xVt induced by the partition V of the state space: x^qx' if and only 
if flow (x) = flow{x'). From Theorem[T]it follows that '^n is proposition-preserving (3. 
in the sense that for all x,x' £ D and for all tt G 11, a; ^ tt if and only if a;' \= n. 

The discrete or qualitative abstraction of a PA transition system S-TS, called qualitative 
PA transition system, is now defined as the quotient transition system of S-TS, given the 
equivalence relation '^o 
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Definition 4 (Qualitative PA transition system). The qualitative PA transition system 
corresponding to the PA transition system S-TS = (il, L, 11, |=) is S-QTS = 

(^/ ~si:L, n, ^^fj : l=~n ) ■ 

Proposition 2 (Qualitative PA transition system). Let S-QTS = (fi/^fj, L, 11, H~q) 
be the qualitative PA transition system corresponding to the PA transition system S-TS = 
(r!,L,n,~.,h). Then 

1. n/^,, = V- 

2. ^^(j C V X L X V, such that D -^.^^-^ D' if and only if there exists ^(i) G Ss and 
T, r', < T < t', such that ^(r) G D, C(t') G D' , and 

• if ; = mi, then for all t G [r, r']: ^(t) e D = D', 

• iil = dim+, then for all t G {t,t']: £,{t) e D' ^ D, 

• iil = dim,-, then for all t G [r, r'): ^(i) e D ^ D'; 

3. |=~n^ P X n, such that Z? |= Dsign = 5* if and only if for all x £ D, S{x) — S. 

Proof. First, the quotient space il/^^ equals V by the definition of the equivalence relation 
The second part of the proposition follows from the definition of as a relation 

on ^/r^a X L X il/^s^ such that D ^Ls^ D' if and only if there exist x £ D and x' G D' 

such that X x' The expression for — ^^f, is a direct consequence of il/^j-, = 23 and 
Definition [21 Third, \=r^^ is a relation defined on i^/r^^i x IT such that D \=,^^ it if and only 
if there exists x £ D such that a; ^ tt Q. The properties considered here are of the type 
Dsign = S (Definition [3]) . Theorem [T] implies the invariance of S for all a; G Z?. □ 

Notice that the transitions labeled by dim^ or dini^ connect two different flow domains, 
since in Proposition [2] we require that D ^ D' . This corresponds to a continuous evolution 
of the system along which it switches from one flow domain to another. On the contrary, the 
transitions labeled by int correspond to the continuous evolution of the system in a single 
flow domain. 

As for E-TS, we define any sequence of fiow domains (I?°, . . . , Z?™), m > 0, as a 
run of S-QTS if and only if for all j G {0,...,m — 1}, there exists I G L such that 

—^-^si 13^+^. The satisfaction relation ^^f, associates to every run a qualitative de- 
scription of the evolution of the derivatives over time. E-QTS can be represented by a 
directed graph G = (P, ^^j^), called the state transition graph. The paths in G represent 
the runs of the system. The state transition graph corresponding to the two-gene example 
is represented in Figure [3l[b) , and (I?^'^, D^-^ , D^-^, D'^-^, D^-^) is an example of a run. 

It directly follows from the definitions of quotient transition system and simulation of 
transition systems that S-QTS is a simulation of S-TS. The converse is not true in 

general, because S-QTS and E-TS are not bisimilar. 

Proposition 3. S-QTS is a simulation of S-TS. 
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As a consequence of Proposition [Sj if there exists a run (x", . . . , x™) of S-TS, then there 
also exists a run (£)", . . . , Z?™) of S-QTS such that e D\ for alH G {0, . . . , m}. In other 
words, E-QTS is a conservative approximation of S-TS. 

4.2 Invariance of qualitative PA transition systems in parameter 
space 

S-QTS provides a quahtative picture of the dynamics of a genetic regulatory network. This 
picture generally depends on the values of the parameters in the PADE model. Since exact 
numerical values for the thresholds and the production and degradation constants k and 
7 are usually not available, it is important to verify to which extent S-QTS is robust to 
variations in the parameter values. 

We therefore introduce a second equivalence relation ~r ^ FxT, defined on the parameter 
space r C R?^Q of the PA system, with q the number of parameters in S. Two parameter 
vectors p,p' € F are equivalent, if their corresponding qualitative PA transition systems, and 
hence their state transition graphs, are isomorphic. Given the equivalence relation ~r, we 
denote by F/^p the quotient parameter space. That is, F/^^ is a partition of the parameter 
space consisting of sets over which the qualitative PA transition system is invariant. 

For our purpose, subsets of F defined by the following parameter inequality constraints 
are particularly interesting. 

Definition 5 (Parameter inequality constraints). The parameter inequality constraints of S 
are a set of total strict orders on {0}, . . . , 6*^' } U {ipi{M) \ M £ A4r}, for every « G {1, . . . , rt}. 

The following theorem states that S-QTS is invariant over the subsets of F defined by 
the inequality constraints of Definition [5l That is, the qualitative dynamics of S is robust to 
changes in parameter values that do not change the total order of the thresholds and focal 
concentrations. 

Theorem 2 (Invariance) . Let P C F be a set defined by the parameter inequality constraints 
of S. Then, there exists some Q G F/^p such that P C Q. 

Proof. Let p,p' G P be two parameter vectors. Given that p and p' satisfy the same 
parameter inequality constraints, they lead to the same mode and fiow domain partitions 
(Definitions [1] and [2]) , and thus the same quotient space ^/^^a = T). As will be seen in the 
next section, the relations ^^^^ and can be characterized by means of rules involving 
inequality constraints (Propositions [H to [H) . It can be directly verified that the necessary 
and sufficient conditions in the rules apply in exactly the same way to p and p' . □ 

Suppose that protein A inhibits the expression of gene 6 at a concentration lower than 
that required for the inhibition of the expression of its own gene. This gives the inequality 
constraint 6*^ < 0^ in Figure[3ljc). Moreover, if we further assume that when gene a is active, 
the concentration of protein A tends towards a level above which autoinhibition occurs, we 
obtain the inequality constraint Ka/ja > 0a- For these inequality constraints, and similar 
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constraints for protein B, the state transition graph in Figure [3ljb) is invariant. Whereas 
exact numerical values for the parameters are usually not available, the weaker information 
required for the formulation of the inequality constraints can often be obtained from the 
experimental literature. This will be illustrated in Section [6] for the E. coli example. 

5 Symbolic computation of qualitative PA transition sys- 
tem 

The computation of the qualitative PA transition system S-QTS is greatly simplified by the 
fact that the domains D and the focal sets ^'(Af) are hyperrectangular, which allows them 
to be expressed as product sets, i.e. D = Di x . . . x Dn and ^'(M) = ^'i(M) x . . . x vI/„(M). 
As a consequence, the computation can be carried out for each dimension separately. In 
this section, we will describe rules to determine the set of states V, the satisfaction relation 
) ^-ncl the transition relation — , as well as their implementation in the computer tool 
GNA. In practice the computations reduce to simple checks of ordering relations, which can 
be carried out symbolically by means of the parameter inequality constraints. 

5.1 Computation of states 

In order to compute the states of S-QTS, we need to determine the fiow domain partition 
P of O (Proposition This requires a total ordering of the threshold concentrations 
{9}, . . . ,6'f } and the focal concentrations {ipi{M) \ M G Mr}, i£{l,...,n} (Definition [2]). 
The parameter inequality constraints of S provide this information. 

Each flow domain has associated to it a number of properties, deflned by the satisfaction 
relation \=n..n- From Proposition [2] it follows that computing the relation |=^q amounts to 
computing S{x). The following two propositions, which are direct consequences of the basic 
mathematical properties of solutions of PA systems discussed in Section [Sj show how to 
compute the derivative sign patterns. 

Proposition 4 (Computation of Dsign). Let D E V he an instantaneous flow domain. 

D ^^j, Dsign = {}. 

Proof. There exists no solution remaining in D for some time, so by the deflnition of S{x), 
we have S{x) = {}, for every x € D. □ 

Proposition 5 (Computation of Dsign). Let D G P be a persistent flow domain. 

D Dsign = S ^ {}, S ^ Si x . . . x Sn, and for every i G {1, . . . , n}, 

if Di is located in a threshold or focal concentration hyperplane, then Si = {0}; 
otherwise, 

• — 1 G Si, if and only if for all a; G D there exists tp G '^{mode{D)) such that 
ipi - Xi < 0; 
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• e S'i, if and only if for all a; G D there exists V' G ^{mode{D)) such that 

ipi- 0; 

• 1 e 5i, if and only if for all a; G D there exists S 4'(mode(£>)) such that 
■01 - > 0. 

Proof. Let M = mode{D). Because D is persistent, there exists a solution S^{t) e 
remaining in D for some time. Let ^(r) = a; £ _D for some r > 0. From (fT2l) we infer 
S{x) ^ {}. By Theorem [1] this must hold for all x£D,bo {}. 

Let Di be located in a threshold or focal concentration hyperplane. For every solution 
£,(t) in D, such that ^(r) = a; e D for some r > 0, the derivative of ^(r) exists, and 
({t) e ^^(^(t)), it holds that = and hence S^{x) = {0}. By Theorem [1] this must 

hold for all x G D, so Si = {0}. 

Alternatively, Di is not located in a threshold or focal concentration hyperplane. We only 
consider the case — 1 G S'i (the argument for and 1 is similar). We first prove necessity of 
the inequality condition. —1 £ Si means that for all x G D, there exists a solution £,{t) £ Ss 
with ^(r) — X and t > 0, such that ^^(t) < (Proposition [2] and lfT2|) ). If M is regular, then 
H^i^ir)) = {z^f (V'.(Af) - e»(T))}. Since C^(t) £ ff.(C(^)), it follows that Vi(M) - a;, < 0. 
If M is singular, then H,{^{t)) = T^{{i^^'' {ip^{M') - ^^(r)) | M' £ i?(M)}). ^^(t) < 
impHes that there exists some iplM'), M' £ R{M), such that tpi{M') — Xi <0. Conversely, 
by Lemmas [T] and H] there exist a solution ^ (t) £ Sx; in M monotonically converging from x 
to ip, with ^(0) = x. Because ipi — Xi < 0, we have ^i{t) < 0, t > 0, and thus —1 £ Si{x). 
By Theorem [1] it follows that —1 £ Si, which proves sufficiency. □ 

Notice that the total ordering on the threshold and focal concentrations expressed by the 
parameter inequality constraints (Definition [5]) allows one to decide which of the conditions 
in the second part of Proposition [5] are satisfied. As a prerequisite, '^{mode{D)) needs 
to be determined, which is also straightforward, even for singular mode domains, given 
the definition of the focal set and the parameter inequality constraints. We illustrate the 
application of Propositions |4] and [5] to our two-gene example. 

First, consider the fiow domain D^-^ in Figure IDJa). "^{mode{D^-^)) equals 
{{Ka/ja, i^b/jb)'}, SO that D^'^ IS persistent. We therefore apply Proposition [5] to deter- 
mine S = Sa X Sb, such that D^-^ \=~!:i Dsign = S. Given that Dl;^ = [Q,9\) and 
ipai'm,ode{D^-^)) = Ka/ja, and 6}^ < Ka/ja according to the parameter inequality constraints 
in Figure [3l[c) , it follows that — > 0, for all a; £ Z?^'^ and V G ^'(moc!e(_D^-^)). Conse- 
quently, we infer Sa = {1}. Repeating this procedure in the a;6-dimension, we similarly find 
Sb — {!}, so that D^-^ Dsign — {(1, 1)'}. This means that the solutions in D^-^ are 

strictly increasing in both dimensions. 

As a second example, consider the fiow domain Z3^-^, represented in Figure [ll[b) . = 
mode{D^ '^) is a singular mode domain, so we first have to compute \I>(M^). It can be 
immediately verified that R{M'^) = {M^,M^}, where ^p{M^) = {Ka/la,OY and ^{M^) = 
(0,0)' (Figure EJb)). As a consequence, reci{{il;{M^),ilj{M^)}) = [0,Ka/ja] x {0}. From 
the parameter inequality constraints in Figure [3ljc) it follows that < 6*^ < Ha/ja, so that 
supp{M^) and recI({V'(M3),-(/)(M5)}) intersect at ^(M^) = {(6*2,0)'}. D^ '^ is persistent, so 
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Figure 4: Derivative signs of flow domains: (a) Dsign{D^-^) 
DsigniD"^-^) = {(0, -1)'} and Dsign{D^-^) = {}. 



{(1,1)'} and (b) 



that Proposition[5]applies. We determine S = Sa^Sb, such that D^'^ H~n Dsign = S. Since 
Z?^'^ coincides with a threshold hyperplane, Sa — {0}. Bearing in mind that D'^-^ ~ {0,db) 
and TpbiM'^) = 0, we flnd Sb = {-1}. This results in D"^ "^ Dsign = {(0, -1)'}, which is 
of course consistent with the fact that the solutions sliding along f*-^ monotonicaly converge 
towards ^'(M'*), and are therefore strictly decreasing in the Xfc-dimension. 

5.2 Computation of transitions between states 

In Section |4!T] we have distinguished three types of transitions: int, dim~ , and dim^. We 
will formulate, for each of these three cases, transition rules that can be appHed by means 
of the parameter inequality constraints of S. 

The transitions of type int are easy to determine, since by Proposition [2] they are nec- 
essarily self-transitions, from a flow domain D to itself, which occur if and only if D is 
persistent. 

Proposition 6 (Computation of int transition). Let D ^V. D -^r^„ D if and only if D 
is persistent. 

Proof. By Proposition [2] an int transition occurs if and only if there exist solutions ^{t) G Ss 
remaining in D. That is, if and only if D is persistent. □ 

Recall that the persistence of a domain can be determined by checking whether 
'^{mode{D)) ^ {} (Lemmas [T] and [H) ■ In our two-gene example, the focal set '^{M^) is 
not empty. Therefore, D^-^ is persistent and there exists an int transition on D^-^ . On the 
other hand, \E'(M^) is empty, so D^-^ is instantaneous and does not have an int transition 
(Figure Hb)). 

A dim'^ transition D ''-^ D' is dimension increasing, and therefore requires that 
D C dD' (Section 14. ip . In order to make D' reachable from D, the solutions in D' must 
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point away from D in the dimensions i e {1, . . . , n} for which Di C dD'^. This is expressed 
by the following proposition. 

Proposition 7 (Computation of dim~^ transition). Let D,D' E T> and D C dD'. 

D '^-^^ D' if and only if '^{'mode{D')) ^ {} and there exist x € D, x' E D' , and 
'ip' e ■^{m,ode{D')), such that for alH S {1, . . . , n} for which A C c^D^, it holds that 

(■)/',' - Xi){x[ - Xi) > 0. (13) 

Proof. Let M' = mode{D'). We first prove necessity by contraposition. If 5'(A/') = {}, 
then no solutions remain in D' and a transition from D to D' is not possible. Otherwise, 
suppose that for all x e D, a;' € D', and i/)' G ^'(M'), there exists some i £ {1, . . . , n} for 
which Z3i C dD[ and (■(/'i' — Xi){x[ — Xi) < 0. Furthermore, assume x[ — Xi > for all x € D, 
x' E D' (the case x'.^ ~ Xi < goes analogously). As a consequence, ■i/'i <^ 2;^ < a;^, for all 
X E D,x' £ D', and ?/;' e ^(A/')- By Lemmas [Hand [31 for all solutions iit) E 5s in D' , ^.{t) 
monotonically converges towards ^i{M'). As a consequence, no solution enters D' from D 

and there does not exist a transition D '^-^^ I?'. 

Next, we prove sufficiency. Suppose \I/(M') ^ {} and there exist x E D, x' E D' , and 
ip' E *(M'), such that for alH G {1, . . . , n} for which A C az?-, it holds that {ipl - Xi){x[ - 
Xi) > 0. By the definition of the fiow domain partition (Definition [2]) this holds for all x S D 
and x' E D' . We further assume that x'^ — Xi > for all a; S Z?, a;' ED' (the case x'^—Xi < 
goes analogously). As a consequence, a;i < a;^ < V'i for all a; G -D and x' E D' , and some 
Ip' E *(M'). By Lemma [H for all x' E D' there exists a solution ^(i) E Ss in A, with 
^(0) = x\ which monotonically converges towards ip' . More precisely, by the construction of 
£^{t) in the proof of Lemma [H we have £^i{t) > as long as S,i{t) < As a consequence, it 

is possible for solutions to enter A from D and there exists a transition D D' . □ 

The total strict ordering defined by the parameter inequality constraints allows the nec- 
essary and sufficient conditions for the existence of a dirn^ transition to be straightforwardly 
tested. 




>I'(mode(D3-2)) 
max„ 



Figure 5: c?im+ transitions between flow domains. Representation of the state space, with 
a trajectory entering D^-^ from D^-^ . 
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As an illustration of the proposition, consider the dim"*" transition from D^-^ to D^-^ 
(Figure [H). Notice that D^ "^ e dD^ '^. Moreover, D^-^ g dDl'^, so that x'^ - Xa > 0, for 
all X e Z)^-^ and x' G I?^-^. Since ■ipa{'mode{D^ '^)) = Ka/ja, and 0^ < Ka/7a due to the 
parameter inequality constraints (Figure [3ljc)), we have — x'^ > 0, for all x' € D^'^ and 

ip' e '^{mode{D^-^)). As a consequence. Proposition [7] allows us to infer that D^-^ '^^^^ r^^^ 
D^-^. In a similar way, we can infer that there is a dim~^ transition from D^-^ to D^-^. 
However, as expected from the direction of the flow in E, a dim~^ transition from D^-^ to 
D^ '^ is excluded, since ip'f, - x'^ < and x'f, - Xf, > 0, for all x G D'^-^ , x' G D^ '^ and 
iP' G -^{modeiD^-^)). 

The rule for dim~ transitions is almost symmetric and given in Appendix [Al 



5.3 Computer implementation 

In summary, given a FADE system E and parameter inequality constraints, we can compute 
its qualitative abstraction, that is, the qualitative PA transition system E-QTS, by means 
of Propositions [4] to m Instead of numerically computing the derivative signs in the domains 
and the transitions between domains, we have developed symbolic algorithms exploiting the 
parameter inequality constraints of Definition [H In particular, we map the total strict order 
on B, = {el,...,0P'}U{^,{M) I M G Mr}, i G {1,. to the set = {0, . . . , + 1}. 

This allows Di and 5'i(M) to be expressed as intervals on C. The conditions in the proposi- 
tions then naturally translate into simple inequality tests on integer coordinates of domains 
and focal concentrations (see 0| for details). 

The computation of E-QTS has been implemented in a new version of the computer tool 
Genetic Network Analyzer (GNA) In order to facilitate the analysis of E-QTS, the state 
transition graph generated by GNA can be exported to standard model-checking tools Hke 
NuSMV and CADF 

In practice, since the number of flow domains in the state space grows exponentially 
with the number of genes in the network, it is not usually possible to compute the complete 
state transition graph. However, knowledge of the part of the graph reachable from a (set 
of) initial flow domain(s) is often sufficient to answer most of the questions of biological 
interest. In GNA it is possible to either compute the complete qualitative FA transition 
system or carry out a reachability analysis from a specified set of initial domains. Moreover, 
GNA identifies the equilibrium states, that is, the fiow domains corresponding to (sets of) 
equilibrium points, by testing whether D \=^^-i G Dsign. For each of the equilibrium states, 
the attractor set is computed, that is, the set of states from which the equilibrium state is 
reachable. 
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6 Qualitative analysis of nutritional stress response in E. 
coll 



In case of nutritional stress, an Escherichia coli population abandons exponential growth 
and enters a non-growth state called stationary phase. This growth-phase transition is 
accompanied by numerous physiological changes in the bacteria, concerning among other 
things the morphology and the metabolism of the cell, as well as gene expression |37|. On 
the molecular level, the transition from exponential phase to stationary phase is controlled 
by a complex genetic regulatory network integrating various environmental signals. The 
molecular basis of the adaptation of the growth of E. coli to nutritional stress conditions 
has been the focus of extensive studies for decades [s^l- However, notwithstanding the large 
amount of information accumulated on the genes, proteins, and other molecules known to be 
involved in the stress adaptation process, there is currently no global understanding of how 
the response of the cell emerges from the network of molecular interactions. This suggest 
the use of modehng and simulation tools to study the dynamics of the stress response. 
However, with some exceptions numerical values for the parameters characterizing the 
interactions and the molecular concentrations are absent, which makes it difficult to apply 
traditional methods. 
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Figure 6: Network of key genes, proteins, and regulatory interactions involved in the nutri- 
tional stress network in Escherichia coli. The contents of the boxes labelled 'activation' and 
'supercoiling' are detailed in [4j|. 



The above circumstances have motivated the qualitative analysis of the nutritional stress 
response network in E. coli by means of the method presented in this paper fi^]. On the 
basis of literature data, we have decided to focus, as a first step, on a network of six genes 
that are believed to play a key role in the carbon starvation response (Figure [6]) . The 
network includes genes involved in the transduction of the nutritional stress signal (the 
global regulator crp and the adenylate cyclase cya), metabolism (the global regulator fis), 
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cellular growth (the rrn genes coding for stable RNAs), and DNA supercoiling, an important 
modulator of gene expression (the topoisomerase top A and the gyrase gyrAB). 

Based on this information, a FADE model of seven variables has been constructed, one 
protein concentration variable for each of the six genes and one input variable {u signal) 
representing the presence or absence of a carbon starvation signal [4J|. As an illustration, 
the piecewise-afHne differential equation and the parameter inequality constraints for the 
state variable XtopA are given below. 

itopA — l^topA S {x gyrAB J d gyrAB top A, OtopA 'J top A XtopA 

< (^topA < (^lopA < l^topA/ltopA < max top A 

The equation and inequalities state that in the presence of a high concentration of Fis 
{s^{xfis,Ojj^^) = 1), and of a high level of DNA supercoiling {s^ {x gyrAB, S'j^yrAB) 
■ s~(xtop>i, ^top^) = 1), the concentration of Top A increases, converging towards a high 

value {ntopA/ltopA > OIopa)- 

Using the computer tool GNA, we have performed reachability analyses on the qualitative 
FA transition system associated with the FADE model. The simulation of the entry into 
stationary phase has given rise to a state transition graph of 66 states, computed in less 
than 1 s on a FC (800MHz, 256Mb). Figure [7] represents the temporal evolution of two of 
the protein concentrations in a run. The predicted expression profiles are consistent with 
the observations [43 |. 

The coupling of GNA with model-checking tools [sj has allowed a more systematic veri- 
fication of observed dynamical properties. In the measured concentration of the global 
regulator Fis is shown to decrease and become steady in stationary phase, which is charac- 
terized by a low concentration Xrm of stable RNAs. This can be expressed by means of the 
following CTL formula: 

'EF{±fis < A EF(i/is = A Xrm < drm))- (14) 

The verification of lfT4|) takes a fraction of a second to complete and shows that the for- 
mula holds for the state transition graph. The observation that the transcription of cya is 
negatively regulated by cAMP and CRP [s^l is also reproduced by the model. In fact, the 
following CTL formula is satisfied by the graph: 

AG(Xcrp ^ ^ crp ^ ^cya ^ ^ cya ^ ^signal ^ ^signal ^ EF Xcya ^ 0)- (f^) 

The formula expresses that always, for high levels of CRF and Cya, in the presence of the 
carbon starvation signal, the system can eventually reach a state in which the expression of 
cya decreases. 

The application of the method has led to new insights into how the carbon starvation 
signal results in the slowing-down of bacterial growth characteristic for the stationary phase 



44| . In summary, the analysis has brought to the fore the role of the mutual inhibition of 



Fis and CRP, which in the presence of a carbon starvation signal results in the inhibition of 
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Figure 7: Temporal evolution of Fis and CRP concentrations in the run (D^, . . . , D^"^). For 
every in the run, the profiles display the corresponding concentration intervals Dj-^^ and 
Dirp- The symbols t, i, and o indicate the sign of the derivative for persistent states. 



fis and in the activation of crp. This causes a decrease of the expression of the rrn genes, 
which code for stable RNAs and are a rehable indicator of cellular growth. In addition 
to this increased understanding of the transition from exponential to stationary phase, the 
model has yielded predictions on the occurrence of damped oscillations in some of the protein 
concentrations after a nutrient upshift, predictions that are being tested in our laboratory. 

We are currently working on extended models of the nutritional stress response network. 
The increase in the number of variables naturally leads to the generation of larger state 
transition graphs by our method. In order to investigate the upscalability of the method 
more systematically, we have applied it to a FADE model with nine state and two input 
variables, describing the initiation of sporulation in the bacterium Bacillus subtilis ji^l . More 
specifically, we have analyzed the response of the cell to carbon starvation in the case of the 
wild-type and a dozen of mutant strains. On average, the state transition graphs generated 
by our method consist of 585 states, with a maximum of 2234 states (computed in 14 s). 
The analysis of graphs of this size does not pose any problems for current model-checking 
tools, which shows that our approach is upscalable to large and complex networks. 
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7 Discussion and conclusions 

We have presented a method for the qualitative analysis of hybrid models of genetic regula- 
tory networks. The method is based on a class of piecewise-affine differential equation models 
that has been well-studied in mathematical biology. By defining a qualitative abstraction 
preserving the sign pattern of the derivatives of concentration variables, the continuous PA 
transition system associated with a PADE model is transformed into a discrete or qualitative 
PA transition system whose properties can be analyzed by means of classical model-checking 
tools. The qualitative PA transition system is a conservative approximation of the under- 
lying continuous PA transition system and can be easily computed in a symbolic manner 
by exploiting inequality constraints on the parameters. We have applied the implementa- 
tion of the method to the analysis of a system whose functioning is not well-understood by 
biologists today, the nutritional stress response in the bacterium E. coli. 
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Figure 8: (a) Two-dimensional projection of a sHce of the phase space of the E. coli stress re- 
sponse model for the variables X and Xcya' (b)-(c) Partitioning into (b) mode domains and 
(c) flow domains of the projection, (d)-(e) Excerpts of the state transition graph resulting 
from the qualitative abstraction based on (d) mode domains and (e) flow domains. 

The results of this paper extend our previous work on the qualitative analysis of PADE 
models of genetic regulatory networks [2l|, [i^l- In particular, we have defined a refined 
partitioning of the state space which underHes a qualitative abstraction preserving stronger 
properties of the qualitative dynamics of the system, i.e. the derivative sign pattern. The 
resulting qualitative PA transition system is better adapted to the abstraction level of the 
experimental data, in the sense that it avoids verification of dynamical properties to be 
over-conservative. Consider Figure [HI which compares two-dimensional projections of a 
state-space sHce of the stress response model. Depending on whether mode domains or fiow 
domains are used as the abstraction criterion, the state transition graph will be different 
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(compare (d) and (e) of Figure[8]). Whereas the CTL formula EF{xcrp > AEF {xcrp < 0)) 
holds for the graph in (d), this is not true in (e), thus reveahng that the coarse-grained 
abstraction may cause models to escape refutation by available experimental data. Judging 
from our experience with several FADE models of bacterial regulatory networks, the use of 
a fine- grai ned abstractions leads to only a modest increase in the size of the state transition 
graph [10|. That is, the increase in precision does not exclude the application of the refined 
abstraction to larger systems. 

The hybrid character of the dynamics of genetic regulatory networks has stimulated the 
interest in the appHcation of hybrid-systems methods and tools over the past few years [l|, [TH, 



21|, |2a, [301 . Our approach differs from this related work on several counts. Whereas we use 



piecewise-affine deterministic models, other groups have employed multi-affine and related 



deterministic models [Uj, |42| or stochastic models [3a|. These models are less restrictive 
and thus provide a more precise description of the network interactions. However, they 
are more difficult to analyze and in higher dimensions usually require the appHcation of 
numerical techniques. This is not straightforward to achieve for most biological systems, 
since numerical information on parameter values is usually imprecise or simply not available. 

The FADE model s ifll) i n this paper have been well-studied in mathematical biology 
3 H m, S 0, [iaH, M, \M, SiS, S, SI , and have also formed the basis for other 



work in the field of hybrid systems [28|. However, contrary to [23|, we take into account 
the dynamics of the system on threshold hyperplanes, where equihbrium points and other 
phenomena of interest may occur [il,[33 (but see [H for ideas on how to extend the approach 
in [i^l). In [2^ the partition of the state space is dynamically refined, so as to arrive at a 
qualitative PA transition system that is a better approximation of the original PA system. 
This requires the use of quantifier elimination methods jsH] which allow to decide more 
general and more powerful properties than the rules proposed in Section [5] of this paper, but 
that also incur higher computational costs. The approach of this paper allows us to fully 
exploit the favorable mathematical properties of the FADE models llj , and thus promote 
the upscalability of the method to large and complex networks (Section [6]) . 

From a more general perspective, our approach can be seen as an application of the 
notion of discrete abstraction, commonly used to study the dynamics of systems with an 
infinite number of states 0, [1, 0, 0, S| • Much work has focused on the identification 
of classes of continuous and hybrid dynamical systems for which bisimulation relations with 
finite transition systems are guaranteed to exist. The results of this paper can be seen as 
showing that the weaker simulation relation may also be of considerable practical interest, 
especially for classes of systems for which the existence of a finite bisimulation cannot be 
guaranteed. Discrete abstraction criteria similar to the one used in this paper, based on 
the sign of the (higher) derivatives of continuous variables, have been proposed by other 
authors in the fields of hybrids systems 0, IZ^ and control theory and qualitative reasoning 
In comparison with these approaches, our work deals with a less general class 



of models. However, this allows the development and implementation of efficient, tailored 
algorithms for the practical computation of the qualitative dynamics of the system, even on 
(intersections of) threshold hyperplanes, where discontinuities may occur. 
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The possibility to use efficient algorithms for the computation of the qualitative PA 
transition system rests, to a large extent, on the approximation of the set K{x) in (0]) 
by the set H{x) in Because the latter set is hyperrectangular, the computation of 

domains, transitions, and sign patterns can be carried out seperately in every dimension, 
using the ordering of parameter values fixed by inequality constraints. Because H{x) is an 
overapproximation of K{x), the state transition graph may contain sequences of states that 
would not occur in the graph obtained by using K{x). As a consequence, a FADE model 
may fail to be rejected by an observed time-series of measurements of the concentration 
variables. However, due to the fact that the approximation of H{x) by K{x) is conservative, 
a FADE model will never be falsely rejected. An obvious direction for further research would 
be to see whether finer overapproximations of K{x) can be found that still allow tailored 
symboHc algorithms to be used that do not compromise the upscalability of the method. 
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A Proofs 

Proposition 1 (Reachability equivalence). For all x, x' G £7, there exists a solution ^(t) £ S 
and T, t', such that < r < r', ^(r) = a;, and ^(r') = x' if and only if there exists a run 
(a;°, . . . , X™) of S-TS such that x° = x and = x'. 

Proof. We first prove necessity. Consider two points x,x' e fi, a solution ^(i) e Ss, and 
T, r' such that < r < t', ^(r) = x, and ^(r') = x' . If r = r', then (C(r)) is a trivial 
run satisfying the conditions of the proposition. Otherwise, t < t' and we denote by 
Z?°, . . . , the time-ordered sequence of fiow domains traversed by £,{t) on the time interval 
[r, t']. £)°,...,D'" is a finite sequence, since by definition any solution of S reaches and 
leaves finitely-many times a threshold hyperplane during a time interval. If m = 0, then by 
Definition [3] there exists an int transition from ^(r) to C(t'), and consequently, {£.{t),^{t')) 
is a run satisfying the conditions of the proposition. Otherwise, m > 0, and we denote by 
(T°, . . . , (T™~^, the switching times, formally defined as cr^ = sup{i G [t, r'] | £^{t) G D^}, 
j e {0, . . . , m — 1}. Finally, we introduce a sequence of time instants t°, . . . , r™ such that 
^{t^) e , j e {0,...,m}. More precisely, we define as t, as (cr-'"^ -I- cr^)/'2, for 
all j e {1, ... ,m — 1}, and r™ as r'. It is not difficult to show by induction on j that 
(^(r°), . . . ,^{t^)) is a run, for every j £ {0, . . . , m} 0. 

Next, we prove sufficiency. Consider a run (x'', . . . , x™) of S-TS, with x'^ = x and 
X™ = x' . If m = 0, then x = x' and any solution ^(i) £ Ss with ^(0) = x satisfies the 
conditions of the proposition for t = t' = 0. In the sequel, we suppose to > 0. Then, by 
Definition m there exists a sequence of solutions {S,^{t), . . . ^ind of time instants 
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(r", . . . , r™), such that for all j G {0, . . . m— 1}, £,^t) is defined on the time interval [t^ , t^^^], 
with < T^~^^ , and satisfies ^■' (r-' ) = and £,^{t^^^) — x^^^. It can be straightforwardly 
shown j3] that the concatenation of the solutions i^{t), j G {0, . . . , m — 1}, is a solution that 
satisfies the conditions of the proposition. □ 

Proposition 8 (Computation of dim~ transition). Let D,D' E V and D' C dD. 

D '^^'^ D' if and only if "^[mode{D)) ^ {} and there exist x E D, x' <E D' , and 
tp G '^{mode{D)), such that either (a) for all i G {1, . . . ,n} for which D^ C dDi, it holds 
that 

{i^i - x'i){x'^ - X^) > 0, (16) 

or (b) it holds that ip = x' . 

Proof. Let M = mode{D). We first prove necessity by contraposition. If ^'(M) = {}, then 
there exists no solution remaining in D for some time, and a transition from D to D' is 
not possible. Otherwise, suppose that for all a; G D, x' G D\ and V' G ^(-W), there exists 
some i G {1, . . . ,n} for which C dDi, [ipi — x[){x[ — Xi) < 0, and tp ^ x' . We assume 
X ■ — > for all X E D, x' e D' (the case a; ■ — < goes analogously). 

If the inequality is strict, then ipi < x[, for all x' G D' and ip S ^'(M). By Lemma[3l for all 
solutions ^(i) G Ss in 13, ^^(i) monotonically converges towards '^i(M). As a consequence, 

no solution enters D' from D, and there does not exist a transition D ''-^ If the 

inequality is not strict, then by definition of the fiow domain partition (Definition [2]) , we 
have V'i = x[ for all x' G D' and "0 G \1/(M). From the same definition it also follows that 
either ipi = ipi{M) (if AI is regular) or ipi = maxjv/'gi?(M) ''PiiM') (if M is singular). It can 
be shown by construction in the way of the proof of Lemma H] that solutions reach D' only 
asymptotically, as t oo. Moreover, from the definition of ^'(M) and Lemmas [T] and [31 
it follows for every j G {1, . . . ,n} that either ^j{t) G 5'j(M) or that £,j{t) monotonically 
converges towards \l/j(M), as long as ^{t) remains in D (and thus in M). Therefore, because 
^{t) reaches x' G D' asymptotically, we have x' G ^'(M). This contradicts the assumption 
^ x' , so there does not exist a dim~ transition. 
Next, we prove sufficiency. Suppose 'I'(M) ^ {} and there exist x E D, x' e D', 
and tp G ^'(Af), such that (a) for all i G {l,...,n} for which Di C 913 it holds that 
{ipi — Xi){x[ — Xi) > 0, or (b) it holds that ip — x' . Consider (a) and assume x[ — Xi > Q for 
all X £ D, x' € D' (the case x^ — Xi < goes analogously). As a consequence, Xi < x'i < ipi 
for all a: G -D and x' G D' , and some ip G 5'(Af). By Lemma HJ for all x G £> there 
exists a solution S^{t) G Ss in D, with ^(0) — x, which monotonically converges towards 
tp. We choose x such that the corresponding solution enters D' from D, and thus obtain a 

transition D D' . In case (b) we find Xi < x[ = ipi, and by the same argument there 

exists a solution that enters D' from D. However, in this case the solution only reaches D' 
asymptotically, as i ^ oo (see the proof of necessity). □ 
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